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* Abstract  -  This  paper  describes  a  new  approach  to  numerical  solution  of  the  high  frequency  approximation  to  the  wave  equation. 
Traditional  solutions  to  the  Eikonal  equation  in  high  frequency  acoustics  are  obtained  via  ray  tracing.  In  ray  tracing,  the  numerical 
grid  on  which  solutions  are  computed  becomes  distorted  over  time  as  rays  diverge  from  the  initial  wavefront,  reducing  the  accuracy  of 
the  solution.  The  level  set  method  is  a  fixed  grid  method,  whereby  the  user  controls  the  underlying  grid,  and  hence  the  accuracy,  of  the 
solution. 


I.  Introduction 

Accurate  and  computationally  efficient  simulation  of  physical  processes  is  key  to  reducing  the  need  for  expensive  and 
environmentally  risky  at-sea  system  level  experiments  and  data  collection.  Level  Set  Methods  (LSM)  are  generic,  computational 
techniques  developed  by  Osher  and  Sethian  [1]  for  tracking  the  evolution  of  moving  curves  and  surfaces.  Traditionally, 
propagation  models  use  ray  tracing  to  solve  via  method  of  characteristics.  When  rays  (characteristics)  diverge,  eventually  they  do 
not  cover  enough  of  the  physical  space  so  that  accurate,  well-resolved  solutions  are  not  available  on  any  uniform  grid.  For  this 
reason,  there  has  been  interest  in  developing  methods  for  solving  general  propagation  problems  in  a  fixed  frame  of  reference.  The 
advantage  of  this  approach  is  that  standard  partial  differential  equation  (PDE)  solvers,  e.g.,  finite  difference  approximations,  can 
be  employed  to  solve  the  problem  on  uniform  grids  which  may  then  be  refined  to  produce  higher  accuracy.  LSM  achieve  this  by 
representing  the  propagating  surface  as  the  zero  level  set  of  a  higher  dimensional  function.  This  function,  say  0(x,  t )  ,  is 
designed  so  that  0(x,O)  =0  for  x  on  the  initial  propagating  surface  in  phase  space.  This  zero  level  set  is  then  transported  via 
the  underlying  velocity  field.  In  the  case  of  acoustic  propagation  in  isotropic  media,  the  propagation  direction  is  normal  to  the 
propagating  surface  (wavefront).  LSM  are  attractive  due  to  their  robustness  -  they  easily  handle  the  evolving  topology  of  the 
surface  being  tracked,  the  normal  vector  and  curvature  can  be  extracted  at  any  point  of  the  front  from  the  level  set  function 
(provided  the  normal  and  curvature  are  well-defined  at  that  point),  and  it  is  straightforward  to  extend  the  theory  to  higher 
dimensions. 

Software  packages  based  on  several  computational  approaches  in  addition  to  ray  tracing  already  exist  which  can  accurately 
solve  the  equations  of  acoustic  propagation.  For  instance,  to  compute  the  full  wave  equation  solution,  one  can  use  normal  modes, 
or  Parabolic  Equation  (PE)  methods.  However,  when  computational  expense  is  at  a  premium,  such  as  for  real-time  simulations, 
these  are  not  appropriate  methods  especially  at  high  frequencies  where  required  grid  sizes  become  very  large  (in  one  dimension, 
we  need  at  least  2  points  per  wavelength  to  resolve  the  wave).  Ray  tracing  is  therefore  the  current  standard  for  high  frequency 
propagation  modeling.  LSM  may  provide  an  alternative  to  ray  tracing  for  solving  the  high  frequency  approximation  to  the  wave 
equation  that  allows  the  simulation  user  greater  control  over  the  accuracy  of  the  solutions. 

This  work  builds  upon  the  foundation  established  by  Osher,  Cheng,  Kang,  Shim,  and  Tsai  in  [2]  in  which  a  basic  level  set 
method  for  geometric  optics  is  introduced.  Eulerian  geometric  optics  has  been  a  topic  of  intense  research  in  the  scientific 
computing  community  for  quite  some  time.  Benamou  [3]  provides  an  overview  of  approaches  to  this  problem.  The  most  similar 
to  the  level  set  method  is  the  segment  projection  method  [4]  in  which  wavefronts  are  tracked  in  phase  space  as  projections  onto 
each  two  dimensional  subspace  of  the  three  dimensional  phase  space.  This  method  is  effective,  but  requires  complicated 
bookkeeping  in  order  to  reconstruct  wavefronts.  The  approach  of  Osher,  et  al.  propagates  the  wavefront  in  the  full  phase  space. 
Reference  [5]  builds  upon  [2]  by  extending  the  method  to  propagation  in  anisotropic  materials.  In  [6],  an  efficient  method  for 
incorporating  reflecting  boundaries  is  introduced. 
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In  this  paper,  we  apply  these  foundations  to  the  specific  problem  of  high  frequency  acoustics.  In  section  II,  we  provide  an 
overview  of  the  method  and  a  description  of  the  implementation.  In  section  III,  we  present  some  preliminary  results 
demonstrating  the  algorithm’ s  performance  in  a  few  sample  cases  including  varying  sound  speed  profiles  and  reflecting 
boundaries. 


II.  Description  of  approach 


A.  Level  Set  Methods 

LSM  were  originally  developed  for  solving  Hamilton-Jacobi  type  equations.  A  Hamilton-Jacobi  equation  is  a  first-order 
nonlinear  PDE  having  the  form  ut  (x,  t )  +  H (x,  Vu)  =  0 ,  where  H  is  a  nonlinear  function.  In  the  case  of  the  Eikonal  equation 

derived  from  the  high  frequency  approximation  to  the  wave  equation,  H (x,  Vu)  =  ±c(x)|  V u\  ,  with  the  nonnegative  function 
c(x)  specifying  the  propagation  speed  of  the  medium;  the  sign  ambiguity  refers  to  the  direction  of  propagation  (inward  or  outward 
from  the  initial  condition).  For  simplicity,  we  only  consider  the  case  //(x,  Vu)  =  +c(x)|Vw| . 

Two  difficulties  must  be  addressed  when  solving  the  Eikonal  equation  in  a  fixed  frame  of  reference:  first,  the  solutions  may  be 
multi-valued  (e.g.,  caustics),  and  second,  solving  over  the  entire  spatial  grid  implies  substantially  increased  requirements  for 
storage  and  processing.  The  first  issue  is  resolved  by  considering  the  wavefront  in  a  higher  dimensional  reduced  phase  space.  In 
two  spatial  dimensions,  for  instance,  the  phase  space  is  four-dimensional,  but  this  can  be  reduced  to  three  dimensions  since  only 
the  direction  of  propagation  is  important.  This  means  that  the  wavefront  is  represented  in  space,  time,  and  a  local  phase  direction. 
The  wavefront  is  embedded  in  two  level  set  functions,  cp  and  \j/,  which  evolve  in  phase  space  according  to  a  first  order  system  of 
transport  equations 

\<bt  +V  *V^  =  0 

\  ,  (l) 

{y/t+V -Vy/  =  0 


where  the  velocity  field  for  the  motion  is  given  by 

xr  r  ^  ^  dc  .  -  dc  „l7 

V  =  [c  cos#,  c  sin#, — sm# - cos#]  ,  (2) 

dx  dz 

as  derived  in  [4],  and  c  =  c(x,z)  is  the  medium  sound  speed.  The  phase  space  variable  #is  the  angle  of  the  local  ray  direction  (for 
an  isotropic  medium).  The  wavefront  is  propagated  in  this  higher  dimensional  space  as  a  system  of  evolving  surfaces,  where 
solutions  are  single- valued  (Fig.  1). 
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Figure  1:  Implicit  representation  of  higher-dimensional  wavefront  -  the  wavefront  is  found  at  the 
intersection  of  the  evolving  zero-level-set  surfaces;  the  horizontal  axes  represent  physical  space  while 
the  vertical  corresponds  to  an  extended  phase  space  dimension. 


The  resulting  physical  wavefront  is  recovered  by  finding  the  intersection  of  the  zero  level  set  surfaces  of  cp  and  \| /  and,  and  then 
projecting  the  resulting  curve  onto  the  horizontal  (physical)  plane  as  in  Fig.  2.  Since  higher  dimensional  spaces  are  required, 
memory  issues  must  be  carefully  addressed.  However,  by  selectively  updating  the  system  only  locally  near  the  wavefront  and  by 
making  use  of  efficient  numerical  methods,  the  burden  on  both  storage  and  processor  time  can  be  greatly  reduced.  This  is  the 
basis  for  the  local  level  set  method  [7]. 


Intersection  of 
Level  Sets 


Figure  2:  Construction  of  the  wavefront  -  the  physical  wavefront  is  found  by  projecting  the  curve  of 
intersection  of  the  two  zero  level  set  surfaces  (as  can  be  seen  in  Fig.  1)  from  the  higher  dimensional 
representation  onto  the  two-dimensional  spatial  plane.  As  can  be  seen  here,  the  wavefront  is  circular 
since  this  simple  case  uses  an  isospeed  medium. 
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B.  Implementation 


The  advantage  of  fixed-grid  methods  is  that  they  allow  for  ready  application  of  standard  partial  differential  equation  solvers 
which  have  well-established  convergence  properties,  e.g.,  finite  differences  or  finite  elements.  The  level  set  equations  themselves 
comprise  a  decoupled  system  of  simple  first  order  hyperbolic  problems.  For  such  problems,  we  focus  on  upwind  schemes,  which 
take  advantage  of  knowledge  of  the  direction  of  propagation  given  by  the  known  velocity  field  V.  Here,  we  employ  a  semi- 
Lagrangian  method  to  solve  the  transport  equation.  Cheng  uses  a  similar  method  in  [9]  to  implement  his  generic  three- 
dimensional  level  set  code. 

There  are  two  steps  involved  in  our  implementation.  In  the  first  step,  a  standard  ordinary  differential  equation  (ODE)  solver, 
such  as  a  Runge-Kutta  method,  is  used  to  step  backward  along  the  characteristics,  which  flow  according  to  the  known  velocity 
field  V  given  in  (2),  one  time  step  at  each  grid  point  in  phase  space.  This  gives  an  approximate  solution  to  where  the  value  on  the 
grid  at  the  current  time  came  from.  Since  these  points  are  not  necessarily  in  the  grid  itself,  an  interpolation  method  is  used  to  find 
the  value  of  the  function  at  that  point,  which  provides  the  solution  at  the  current  time  on  the  grid.  For  example  in  our  case,  using 
backward  Euler  for  the  time  stepping,  assuming  a  constant  sound  speed  c  and  time  step  At ,  we  solve  for 

0*  “  <P(xi  ,Zj,Ok,tn)  on  the  grid  ({x,  ,  {z;  ,  {0k  ]  where  x  and  z  are  the  physical  space  variables,  with  <9  as  the 

phase  space  variable,  as  follows: 

x*  =  xt  -  cAt  cos(Ok ) 

Stepl:  Compute  /  z*  =  z-  -cAt  sin(0k) 

V  =  ek 

Step  2:  Set  =  f~\x\z,0*) 

In  step  2,  the  quantity  (f)"  '  (x  ,z  ,  0  )  is  estimated  via  linear  interpolation  on  the  data  j^-1 }.  These  steps  are  applied  to 

both  level  set  functions  (f)  and  If/ .  Further  interpolation  is  then  used  to  find  the  intersection  of  the  zero  level  sets  and  recover  the 
wavefront.  Higher  order  methods  can  be  applied  here,  however  the  stability  result  may  not  hold.  Note  also  that  if  the  speed 
function  c(x,z)  is  nonlinear,  an  iterative  method  such  as  Newton-Raphson  is  required  to  solve  step  1. 

This  approach  has  two  important  benefits.  One  is  that  it  effectively  decouples  the  grid  points  from  one  another  in  the  time 
stepping  portion,  making  this  part  of  the  code  readily  parallelizable,  as  well  as  allowing  great  flexibility  in  the  choice  of  points  at 
which  to  solve.  This  becomes  critical  for  a  three  dimensional  application  where  reduced  phase  space  has  five  dimensions  and 
computation  must  be  focused  near  the  wavefront.  The  other  significant  advantage  is  that  this  method  is  unconditionally  stable  - 
although  this  result  depends  on  the  interpolation  scheme  used  in  step  2.  That  is,  it  does  not  suffer  from  a  time-step  restriction 
(i.e.,  the  Courant-Friedrichs-Levy  stability  condition)  as  do  standard  Eulerian  finite  difference  methods,  and  thus  larger  time  steps 
can  be  taken  without  sacrificing  stability;  this  could  be  critical  for  a  real-time  implementation  as  accuracy  will  have  to  be 
sacrificed  for  speed. 

For  boundary  conditions,  if  using  finite  difference  methods,  one-sided  formulas  can  be  derived  at  the  boundaries  to  the  desired 
order  of  accuracy.  For  the  semi-Lagrangian  method  however,  an  inflow  condition  must  be  defined  for  the  case  when  information 
is  transported  into  the  domain  from  outside.  If  the  boundary  is  purely  reflecting,  all  information  is  within  the  domain.  We 
currently  use  this  type  of  condition  for  the  bottom  and  surface  boundaries;  since  it  is  a  high  frequency  algorithm,  transmission  into 

the  bottom  is  not  of  significant  concern.  The  reflection  condition  is  implemented  by  setting  (j)n  (x.  zb ,  dreb  )  =  (/)n  (x? ,  zb ,  0inc  )  on 

the  boundary  points,  Zb  ,  where  0inc  is  the  incident  angle  computed  via  Snell’s  law  from  the  reflected  angles,  0refl ,  in  the  grid.  It 

is  not  physical  to  assume  reflection  on  the  left  and  right  domain  boundaries,  so  in  order  to  compute  a  solution  on  the  boundary  at 
incoming  angles,  we  find  the  time  and  location  at  which  the  boundary  was  crossed,  and  assume  the  solution  is  constant  beyond 
that  point.  This  is  a  numerical  scheme  that  works  as  long  as  all  of  the  important  physical  action  is  sufficiently  far  away  from  that 
boundary.  One  nice  feature  of  the  level  set  method  is  that  what  happens  to  the  level  set  functions  away  from  the  zero  level  sets 
themselves  is  largely  irrelevant.  As  long  as  the  source  is  inside  the  domain,  approximations  made  at  inflow  do  not  have  a 
significant  effect. 
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III.  Results 


We  next  present  some  preliminary  results  obtained  with  our  algorithm.  All  of  the  plots  were  generated  from  our  Matlab® 
prototype  code,  using  the  semi-Lagrangian  method  to  solve  the  relevant  differential  equations.  The  first  case  (Fig.  3)  uses  an 
isospeed  profile  using  a  sound  speed  c  =  1500  m/s  and  no  boundaries.  The  circular  wavefront  gradually  expands  and  exits  the 
domain.  Fig.  3a  displays  snapshots  of  the  wavefronts  at  a  few  time  steps  and  Fig.  3b  shows  the  rays  extracted  from  the  level  set 
results.  As  expected  in  an  isospeed  medium,  the  rays  appear  as  straight  lines. 
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Figure  3:  a)  Wavefronts  in  an  isospeed  medium  with  no  boundaries;  b)  Rays  extracted 

from  a) 


In  Fig.  4a,  the  wavefronts  are  depicted  for  a  similar  case,  but  this  time  with  a  linear  sound  speed  profile  that  gradually 
increases  from  1500  m/s  to  2000  m/s  (see  Fig.  5).  The  propagation  effects  are  subtle,  but  we  observe  the  wavefront 
expanding  faster  at  greater  depths  thus  distorting  the  wavefront  from  its  original  circular  form.  The  ray  trace  for  this  case, 
Fig.  4b,  shows  the  rays  curving  slightly  downwards  in  the  direction  of  increasing  propagation  speed.  A  few  “outlier”  points 
appear  on  the  ray  trace;  we  suspect  that  these  transient  errors  are  due  to  inaccuracies  in  the  contour  finding  routine  we  are 
using.  They  have  not  affected  the  overall  solution,  so  are  not  indicative  of  any  underlying  instability. 

The  next  examples  include  reflecting  boundaries  at  the  surface  and  bottom  of  the  domain  for  a  simplistic  shallow  water 
model.  Again  we  present  cases  with  an  isospeed  medium  (1500  m/s)  and  with  the  linear  profile  depicted  in  Fig.  5.  In  Fig. 

6a,  the  initial  wavefront  expands  and  eventually  reflects  off  of  the  surface  and  bottom  boundaries,  presenting  a  nice 
symmetry.  In  Fig.  6b,  the  ray  trace  shows  the  correct  pattern  with  straight  lines  reflected  back  symmetrically.  The  linear 
profile  example  is  more  interesting;  again  the  subtle  bending  of  the  rays  is  observable  as  is  a  small  distortion  in  the  expanding 
wavefront  pattern  in  Fig.  7. 
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Figure  4:  a)Wavefronts  in  medium  with  linear  sound  speed  profile  and  no  boundaries; 
b)  Corresponding  ray  trace 


Sound  speed  profile 


Figure  5:  Sound  speed  profile  for  the  linear  sound  speed  case  (Figs.  4,7) 


Figure  6:  a)  Wavefronts  in  isospeed  medium  and  reflecting  boundaries  at  the  surface 
and  bottom;  b)  Corresponding  ray  trace 
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Figure  7:  a)  Wavefronts  in  medium  with  linear  sound  speed  profile  and  reflecting 
boundaries  at  the  surface  and  bottom;  b)  Corresponding  ray  trace 


IV.  Conclusion  and  future  plans 

In  this  work,  we  have  outlined  a  novel  approach  to  solving  the  high  frequency  wave  equation  of  acoustic  propagation.  The 
demonstrations  of  simple  propagation  models  in  section  III  above  suggest  that  this  method  holds  promise.  There  still  is  a  need  to 
address  the  errors  that  appear  in  the  non-constant  sound  speed  case;  future  plans  include  incorporating  more  specialized 
interpolation  algorithms  for  identifying  the  level  set  surfaces  and  their  intersections.  This  will  also  be  key  to  obtaining  good 
results  with  more  general  geometry.  Once  this  portion  of  work  is  complete,  a  thorough  validation  against  known  results  will 
ensue. 

The  goal  of  this  project  is  to  build  a  prototype  software  package  that  will  combine  the  basic  framework  of  LSM  with  realistic 
surface  and  bottom  models,  appropriate  loss  models,  and  will  be  optimized  for  computational  speed.  It  is  also  desirable  to  extend 
this  to  allow  for  stochastic  boundary  conditions  (e.g.,  the  ocean  surface),  and  propagation  through  random  media.  Existing 
simulations  involving  underwater  acoustics  assume  a  deterministic  acoustic  response,  and  if  any  variability  is  added,  it  is  often 
post-processed,  hence  not  capturing  its  true  nature.  Since  the  variability  originates  in  the  environment,  this  is  where  it  should  be 
accounted  for. 
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